Analysis of noise-induced bistability in Michaelis Menten single-step enzymatic cycle 

Daniel Remondinip Enrico Giampieri, Armando Bazzani, and Gastone Castellani 
Physics Dept. of Bologna University and INFN Bologna 



o 






o 

• I— I 



> 

O 

^, 

o 



X 



Amos Maritan 

Physics Dept. of Padova University 

(Dated: July 25, 2011) 

In this paper we study noise-induced bistability in a specific circuit with many biological im- 
plications, namely a single-step enzymatic cycle described by Michaelis Menten equations with 
quasi-steady state assumption. We study the system both with a Master Equation formalism, and 
with the Fokker-Planck continuous approximation, characterizing the conditions in which the con- 
tinuous approach is a good approximation of the exact discrete model. An analysis of the stationary 
distribution in both cases shows that bimodality can not occur in such a system. We discuss which 
additional requirements can generate stochastic bimodality, by coupling the system with a chem- 
ical reaction involving enzyme production and turnover. This extended system shows a bistable 
behaviour only in specific parameter windows depending on the number of molecules involved, pro- 
viding hints about which should be a feasible system size in order that such a phenomenon could 
be exploited in real biological systems. 

PACS numbers: 



I. INTRODUCTION 

Many biological phenomena (e.g. memory induction, 
chromatin remodeling, cell-fate determination) are re- 
ceiving a great deal of attention in recent years, due to 
an increasing interest in the description of their com- 
plex behaviour by means of basic biochemical circuitry 
Ul, fi Q- The signal transduction machinery is mainly 
based on enzymatic reactions, whose average kinetics 
can be described within the framework initally proposed 
by Michaelis and Menten (MM). The steady state ap- 
proximation of MM model accounts for the majority of 
known enzymatic reactions, and can be adjusted for the 
description of regulatory properties such as cooperativ- 
ity, allostericity and activation/inhibition [13] • The MM 
equations are still valid at small molecule numbers (as it 
frequently happens in real cells) if the microscopic inter- 
pretation is changed correspondingly 0, 0, Il3 but the 
discrete stochastic aspects become predominant, and a 
deterministic or a stochastic continuous model can not 
describe the system in sufficient detail [J, [ij, [l^ . 

A large class of enzymatic reactions controls the re- 
versible addition and removal of phosphoric groups, 
phosphorylation/dephosphorylation reactions catalyzed 
by kinases and phosphatases respectively. The phos- 
pho/dephosphorylation cycle (PdPC) is thus a post- 
translational substrate modification that is central for 
the regulation of several biological processes [^, [ij . 

How these processes can show a bistable behaviour 
0, [U in the presence of fluctuations Q, reflected 
by a bimodal stationary distribution of protein num- 
ber/concentration, is a crucial question for their mod- 



eling. 

The deterministic version of a single PdPC is not 
bistable in general, but it is hypothesized that external 
noise can trigger such a behavour [l4]. We study this cy- 
cle by a Chemical Master Equation approach, and also (in 
the limit of large molecule number) the related Fokker- 
Planck equation. A closed form for the stationary distri- 
bution of the system is obtained for both approaches: we 
show that the system can not have a bimodal stationary 
distribution due to intrinsic fluctuations only, but an ad- 
ditional external noise obtained by a plausible biological 
mechanism (i.e. the coupling of the system with an en- 
zyme production/ activation reaction) can produce such 
feature. We define analytically the conditions in which 
bimodality occurs, as a function of the reaction param- 
eters (kinetic constants) and system size (number of en- 
zyme and substrate molecules) , and we verify our results 
by numerical simulations with a Gillespie algorithm. 



II. THE MODEL 

The PdPC (also referred to as the futile cycle) is com- 
posed by one phosphorylation and one dephosphorylation 
reaction, catalyzed respectively by enzymes Ei and E2'. 
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The deterministic dynamics of this cycle can be described 
via the MM formalism. Assuming a steady-state approxi- 
mation for both enzymatic reactions A and A^ , we obtain 
the following equations: 
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B. Fokker-Planck approximation 
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Imposing the conservation of the total substrate concen- 
tration, let X be the A molecule concentration, we obtain: 



X = Vm2 



l-x 



Km2 + 1- X 
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(3) 



that can be easily shown to have only one solution inside 
the substrate domain (see |14|). 



A. The CME approach 

Starting from the previous equations, a Chemical Mas- 
ter Equation (CME) approach |16i] is introduced to ac- 
count for intrinsic noise (p„ is the A-molecule distri- 
bution function over the possible states n G [0 : iV], 
D+f{n) = f{n + l)-f{n)): 



Pn 
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(4) 



where 
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N is the total number of the substrate molecules, n is the 
the number of A molecules and the MM constants have 
been accordingly scaled: K'j^j — N-Km and Vlj = N-Vm- 
In the hypothesis of fast relaxation times, the stationary 
solution of eq. (HI describes the statistical properties of 
the reaction in Fig. [T] The stationary distribution p^ is 
derived by imposing Pn{t) = 0; excluding the existence 
of a constant current in the system, we get the condition 



Pn 
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D+liip%n) =ln 
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If we define a potential V{n), such that D+V{n) = 
— ln(f;„/r„_|_i), the stationary solution has the Boltzmann 
form 



CME dynamics can be more easily studied by consid- 
ering a continuous approximation by means of a Fokker- 
Planck (FP) equation, with N is sufficiently large to get 
the natural boundary conditions p*(l) = p''{0) ~ 0. Ex- 
panding in power series the Z)+ operator up to second 
order, from eq. ^ we get: 



dp 



{x, t) = -—C{x)p{x, t) + l—D{x)p{x, t) (6) 



where the drift and diffusion coefficients are defined as 

C{x) = g{x) - r(x), D{x) = [r{x) + g{x)] /N, x = n/N 
{0<x< l),pn^pix)/N, and 



9ix) 
r{x) 



Vm2 ■ (1 - a;) . 

Km2 + 1 - X ' 

Vmi ■ X 
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up to an error 0{N^^). We remark that the diffusion 
coefficient becomes negligible for large N with respect to 
the drift term: the fiuctuations scale as 1/a//V according 
to the law of large numbers. The stationary solution 
p^ix) can be written explicitly: 



fix) 
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dy 



that for the symmetric case Vmi = Vm2, Kmi = Km2 = 
Km reduces to 



p'{x) ^ F [Km + 1 - x) [Km + x) {Km + 2a; - 2x^)^'" 

(7) 
where i^ is a normalizing constant. 

The FP stationary solution is a good approximation 
of CME when the drift coefficient is of the same order 
of the diffusion coefficient (i.e. C{x) ~ 0{1/N), see Ap- 
pendix). Therefore, we can apply the FP approximation 
nearby the critical points x* where C{x) ~ 0. By us- 
ing a Gaussian approximation, in the neighborhood of a 
critical point x^ we get the distribution (C' {x) = dC/dx) 



p'^ix) oc exp 



\C'{x,)\ {x-x,f 
D{x^) 2 



(8) 



In the symmetric case, for generic enzyme concentration, 
we explicitly compute 



Pn 
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'V(n) 



(5) 



where i^ is a normalizing constant. According to ([5|), the 
maxima and minima of the distribution are obtained by 
imposing D^V{n) = 0, extending the n domain to the 
set of real numbers. This leads to gn/rn+i — 1, similar to 
Q and with an unique solution inside the [0 : N] domain. 
This result is also confirmed by Gillespie simulations of 
the dynamical process. 
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(9) 



The variance D(a;»)/|C"(a:;*)| scales as {KmN)~^ and, in 
general, is weakly dependent on the enzyme concentra- 
tion. Therefore, in the continuous limit iV ^ 1, we get a 
finite variance only for Km ^ 1. In Fig. [T]we compare 
the exact CME stationary distribution, the exact solu- 
tion of the FP equation and the approximation of Eq. 
(|5]) in the symmetric case, where the condition C{x) <C 1 
is satisfied. In the FP equation (jS]) two minima may ap- 
pear in the physical domain of x for specific values of Km 
only when TV is sufficiently low, but these critical points 
have no meaning being due to a bad approximation of 
the CME. Thus intrinsic noise cannot induce stochastic 
bifurcation in system ([T]). 




FIG. 1: Solutions for the futile cycle: TV = 100, Vm = 1, 
Km = 0.1, El = E2- Squares: CME stationary solution; line: 
FP exact solution; circles: FP approximate solution (|8]). 



When 7=1 (i.e. Ei — E2) we have the trivial solu- 
tion u = (an unique maximum with n = (TV + l)/2, 
X = 1/2), whereas for 7 — 1 > (resp. < 0) w shifts 
towards —a (resp. a). 

Supposing that enzyme concentration can fiuctuate 
around the average value, given ^ = (7 — I) /Km and 
p(^) the corresponding probability distribution, we have 



p{u)du ^ p{^{u)) 



d^ 
du 



du = p{£,{u)) 
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1 



{a + uf 



[a — uY 
(13) 
Under the hypotheses that £, fluctuates around zero 
and p(^) tends suSiciently fast to zero at the boundaries 
(natural boundary condition), we study the conditions 
for bimodality oip{u). The critical points oi p{u) must 
satisfy 



dp{u) d^i 



du 
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(14) 



If we approximate p(^) with a Gaussian distribution (jus- 
tified for a sufficiently large enzyme molecule number and 
K sufficiently small, see Fig. [5]) so that 



p[u) 



1 



1 



[a + uY [a-uY 
the equation (|T4l) reads 






du 



C. Bimodality induced by enzyme noise 

Relaxing the assumption of fixed enzyme concentra- 
tion, we characterize the effect of enzyme fiuctuations on 
the substrate stationary distribution. In the symmetric 
case ifji/i = K'j^.j2 = K'^^and K'(j^ = K'(j2, the equilib- 
rium points of the average equation ([2]) corresponding to 
maxima of p^ can be calculated explicitly as a function 
of the ratio between enzymes 7: 



7 = 



El 
E2 



N -n + l 



K' 



M 



K' 



If one introduces the variable 



n 

N 



N+1 
2N 



n 

N 



M 



u e 



N -n+l 



1 

N 



(10) 



(11) 



where a ~ 1/2 for TV 3> 1, the condition (|T0)) reads 

a — u Km + a + u 



1 = 



a + u Km + a 



where Km = K'j^jN . Assuming that Km ^ 1, so that 
the critical point is quite sensitive to the enzyme con- 
centration, and performing a perturbative approach over 
Kmi the previous equation can be rewritten as 
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= 



If we exclude the symmetric solution u = 0, we get the 
following condition for bimodality (recalling that a^ — 

aj/KM) 



al> 
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9K'^ 



(16) 



The condition (ITBl) can be realized if the enzymatic con- 
centrations fluctuate largely enough around the mean, 
but for the system ([T]) this fluctuation cannot be achieved 
by means of intrinsic noise only. 



III. STOCHASTIC HYPERSENSITIVITY 
INDUCES BISTABILITY 

As we have shown in previous sections, intrinsic noise 
cannot induce bimodality in PdPC, but large enzyme 
fiuctuations can produce this effect for Km sufficiently 
small. This can be achieved by coupling the initial sys- 
tem ([1} with a further reaction involving enzyme activa- 
tion (e.g. by a second messenger like Calcium) or local- 
ization (e.g. inside a delimited region like an organelle). 
Then, in the hypothesis that the activator is in low abun- 
dance (or that the reaction region is small) a situation of 



competition between the two enzymes for the reaction is 
obtained. Under suitable kinetic parameters, the extra 
fluctuation may lead the cyclic reaction alternately to- 
wards one of the two products, producing a bimodal dis- 
tribution function for the substrates. Defining as E^ and 
£"2 the inactive (external to the reaction region) enzyme 
concentrations, and Ea the concentration of activating 
molecules, the full kinetic reaction scheme becomes: 
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These equations are identical to those of the enzymatic 
cycle in ([1} as long as there is a sufficiently large number 
of enzyme molecules, so that the fluctuations in enzyme 
concentration become negligible. A simplifled version of 
the enzyme competition can be obtained by considering 
a direct interchange between the two active enzymes: 

El ^ E2 {E2 ^ Ej" — El j 

with Et total enzyme concentration. The resulting FP 
equation of the system is two-dimensional, and the sta- 
tionary distribution can be expressed as a function of en- 
zyme ratio 7 and substrate concentration x: remember- 
ing that ^ = {j-1)/Km we ha.veps{x,^) = Il{x\(_)ps{S,), 
in which n(a::|^) is the conditional distribution of sub- 
strate given the enzyme concentration. Let us consider 



the case when the PdPC kinetics is much faster than 
the enzyme (achieved by choosing the Kq parameters in 
equations ^ sufficiently large) so that one can apply an 
adiabatic approach to the global system evolution, ap- 
proximating n(x|^) with the stationary solution p^ (x) as 
in ([7]). The substrate distribution Pg{x) is obtained as 
the marginal distribution of ps(a;,^): 



Psix,0 = / Psix\OPsiOd^ = / Ps{x\x^)ps{x^)dx 



(17) 

where we use the univocal relationship between the sta- 
tionary distribution of the enzyme ratio and the equilib- 
rium concentration of the substrate. Applying the calcu- 
lation performed previously, we have shown that Ps(2;») 
is a bimodal distribution if the E1/E2 distribution can 
be well approximated by a Gaussian-like function, and if 
the condition ([T5|) is satisfied. In Fig. [2] we show that 
the Gaussian approximation for the enzyme ratio is valid 
even for a moderate number of enzymes. 

For the simplified model in the symmetric case we get 
a modified version of the previous MM equation (with 
the steady-state approximation and Ei + E2 = 1) 

. Kc (l-E) (l-x) KcEx 

X = 

Km -I- 1 — a: Km + x 

and we have an explicit dependence of the critical point 
x* on enzyme concentration E: 



J 



^{E) = 



2E + Km-1-VI + 2Km -AE + Kjj - 8EKm + 4^^ _^ sKmE^ 

2{2E~ 1) 

I 



(18) 



In Fig. [3] we show the plot of x^,{E) with Km = 0.1 
to point out the sigmoidal behaviour due to the extreme 
sensitivity of the solution to enzyme concentration. Ac- 
cording to (|17p the substrate distribution Pg{x) can be 
computed: 



Psix) oc / exp 
/o 



\c'{x,)\ {x-x,y 

D{x^) 2 



p{x*)dx:, (19) 



where p{x^,) is the critical point distribution defined by 
()13|) and we use the Gaussian approximation for the FP 
stationary distribution. If p(x*) is a bimodal distribu- 
tion (i.e. the condition ((T6| is satisfied), the final dis- 
tribution Pg{x) will also be bimodal when the variance 
D{x^,)/\C'{Xf,)\ is sufficiently small: indeed smaller than 
the distance of the critical points oip{x^). In general, this 
condition is satisfied if {KmN)~^ is sufficiently small, 
which is consistent with the choice A^ ^ 1. In Figure 
m we perform numerical simulations of the CME, using 
a low number of enzyme and substrate molecules. We 
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FIG. 2: Gaussian approximation of p{'y) (7 = E1/E2) for 
El + E2 — 100. Bars: empirical distribution of 7; continuous 
line: gaussian distribution with same mean and variance. 



observe that bimodality occurs only under suitable con- 
ditions depending on the system size: fixing the chem- 
ical reaction constants and the ratio between enzyme 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



FIG. 3: Plot of x,{E) solution ^, Km = 0.1 



FIG. 4: Histograms of a; for TV = 10, 50, 500, obtained by 
numerical simulations run for 5 • 10^ iterations. The ra- 
tio between M (enzyme molecule number) and A'^ (substrate 
molecule number) is set to 0.15; the constants for enzyme pro- 
duction reactions fei, fc_i are set to 2, Kc = 100, Km ~ 0.1. 



is often used instead of the exact Master Equation ap- 
proach, we characterize the conditions in which it ap- 
proximates correctly the system. 

We propose a different approach to generate stochas- 
tic bimodality (that docs not change the number of de- 
terministic stable states, thus a purely noise-driven phe- 
nomenon). With an additional chemical reaction we in- 
troduce competition between enzymes, for binding to an 
activating molecule or for reaching a specific site in which 
reaction can occur (e.g. an organelle), that can be biolog- 
ically feasible. For this system we clarify the conditions 
for which stochastic fluctuations in enzyme concentration 
can lead to bimodality in substrate concentration, and we 
show that it depends on the time scales involved in the 
two reactions (related to the kinetic constants) and also 
on the size of the system (i.e. the number of enzyme and 
substrate molecules involved). 

In particular, if we flx the kinetic constants and the 
ratio between enzyme and substrate molecule numbers, 
bimodality is observed only in specific parameter s, re- 
lated to system size in terms of total number of molecules 
involved. This results define in more detail the feasibil- 
ity of this phenomenon in real biological systems, stating 
that if the biological systems have to exploit noise to 
achieve a bimodal behaviour there must be a relation- 
ship between the chemical parameters of the system and 
its size. 



APPENDIX 



and substrate total molecule number Et/Xt, only in a 
specific interval of molecule number (say Xt, the total 
number of substrate molecules) we have such behaviour. 
The condition (fT6|) can be applied only in the last case 
{N ~ 500) and it is consistent with numerical simulations 
(loss of bimodality) . 

This phenomenon can be explained as follows: in- 
creasing the molecule number, we recover the Michaelis- 
Menten deterministic limit, losing the noise- induced 
bistability; on the other side, if the molecule number is 
too low, even if bistability is not lost from an analytical 
point of view, the separation between the two maxima 
becomes indistinguishable. 



IV. CONCLUSIONS 

In this paper we study the conditions that result in 
a bimodal stationary distribution of a single phosphory- 
lation/dephosphorylation cycle described by Michaelis- 
Menten equations in the quasi-steady state assumption. 
The simple addition of noise (as produced by a Master 
Equation approach) is not enough to achieve bimodal- 
ity in this system that has a unique deterministic stable 
state, and the same result can be stated for the Fokker- 
Planck description of the system obtained as a limit of 
the Master Equation. Since the Fokker-Planck approach 



The FP approximation ([7]) for the stationary solution 
([5]) requires that the condition D+V{n) = —liig„/r„+i 
reduces to 

dV _ C{x) d_ 
dx D{x) dx 

where in the thermodynamics limit N -^ oo where x — 
n/N and r„+i/2 = Nr{x), 5„+i/2 = Ng{x). This is the 
case when (?„ — r„+i is small so that using a perturbative 
expansion we have 
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Therefore if N{r{x) — g{x)) is finite we can approximate 



dV 
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ND+V{n) 



N{r{x)~g{x)) + \^M^)+g{x)) 
r{x) + g{x) 



with an error of order 0{[r{x) — 17(2;)]^). Finally we get 

dV ^^^q{x) — rix) 9 , , , , , ,, 

~ -2N^^ y- + — In r Or + g{x)) 

dx r{x) +g{x) dx ^ ^ ' ^^ " 

Cix) (9 , , , 



The condition Qn — Vn+i «C 1 (i.e. the generation and re- 
combination rates of the enzymatic cycle have almost the 



same value) means that the FP approximation is justified 
only nearby the critical points when gn/^n+i — 1- 
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